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Abstract. In microscopic mechanical systems interactions between elastic structures are often 
mediated by the hydrodynamics of a solvent fluid. At microscopic scales the elastic structures are 
also subject to thermal fluctuations. Stochastic numerical methods are developed based on multigrid 
which allow for the efficient computation of both the hydrodynamic interactions in the presence of 
walls and the thermal fluctuations. The presented stochastic multigrid approach provides efficient 
real-space numerical methods for generating stochastic driving fields with long-range correlations 
consistent with statistical mechanics. The presented approach also allows for the use of spatially 
adaptive meshes in resolving the hydrodynamic interactions. Numerical results are presented which 
show the methods perform in practice with a computational complexity of 0(A'^ log(A'^)). 
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1. Introduction. In many microscopic mechanical systems interactions between elas- 
tic structures are mediated by the hydrodynamics of the surrounding solvent fluid and sub- 
ject to thermal fluctuations. In actual experimental setups and in microscopic devices the 
presence of walls also often plays an important role in mediating the interactions between 
elastic structures. Examples include molecular separation in microfluidic/nanofluidic chan- 
nels [171 118| , manipulation of colloidal probes and oligonucleotides in fluidic assays [6l [9] [19] , 
and the processing of complex fluids in microfluidic devices 26j. In such systems, the hydro- 
dynamic coupling tensors between elastic structures in the flow are no longer homogeneous 
in space and the proximity to the wall plays an important role. For numerical simulations 
this presents a number of challenges. One central challenge is that discretizations no longer 
exhibit translation invariance facilitating a fluid solver based on the Fast Fourier Trans- 
form [101 lll| . When thermal fluctuations are also taken into account, further issues arise 
in generating the required stochastic driving fields consistent with statistical mechanics. 
Stochastic generation methods often rely on factoring the discretized covariance operator of 
the field using the Fast Fourier Transform [3] [2] . For domains represented by spatially adap- 
tive meshes or for domains having a non-rectangular boundary, generation methods based 
on the Fast Fourier Transform can often no longer be used. 

To address these challenges, we introduce new real-space approaches for generating 
the required stochastic driving fields consistent with statistical mechanics. The underlying 
discretization of the hydrodynamic equations is utilized to generate efficiently the stochastic 
driving fields and to account for the hydrodynamic interactions. To efficiently generate 
random variates with long-range correlations, stochastic iterative methods are developed 
which are based on multigrid [8l ll3|[T4ll21| . The stochastic iterative methods introduced for 
fluid-structure systems allow for the efficient generation of stochastic driving fields consistent 
with statistical mechanics both on domains with boundaries and on domains represented 
by spatially adaptive meshes. The ideas presented here are expected to generalize to be 
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Fig. 2.1. Stochastic Eulerian Lagrangian Method (SELM) Coupling. In the SELM approach 
a mixed Eulerian and Lagrangian description is used. The Lagrangian reference frame is used for 
the particles, extended filaments, and elastic bodies. The Eulerian reference frame is used for the 
conservation laws (hydrodynamics). 



applicable more broadly in the development of efficient stochastic numerical methods for the 
simulation of spatially extended stochastic systems. 

2. Stochastic Eulerian- Lagrangian Method for Fluid-Structure Inter- 
actions. To account for the hydrodynamic coupling of elastic structures, we shall use a 
variant of the Stochastic Eulerian-Lagrangian Method (SELM), see [^. Effective equations 
will be derived for the elastic structures and expressed in terms of the differential operators 
of the underlying fluid equations. In the limit in which the fluid is treated as having equili- 
brated to a quasi-steady-state flow, the following effective equations can be derived for the 
elastic structures accounting for their hydrodynamic interactions and thermal fluctuations 

(2.1) ^=i/sELMF + g. 

at 

The F denotes the forces acting on the elastic structure. The g term accounts for the 
thermal fluctuations. The effective hydrodynamic coupling tensor Hselm = Jy"sELM(X) can 
be expressed as 

(2.2) HsELM = -Fpii'^A-^pA. 

It will be assumed throughout that Vx • ^^selm = 0, which ensures phase-space incompress- 
ibility in the dynamics of X. For the stochastic dynamics to exhibit fluctuations consistent 
with statistical mechanics requires the term g for the thermal fluctuations have a covariance 
given bv [21127] 



(2.3) G = (gg^) = 2fesri/sBLM. 

This particular covariance structure can be shown to be a consequence of the principle of 
detailed balance for an ensemble with the Gibbs-Boltzmann distribution when subject to the 



stochastic dynamics given by equation 2.1 see [2]. 

The effective equations [2TJ|23] can be derived as follows. In the quasi-steady-state limit, 
the fluid equations can be expressed as 

(2.4) nAu - Vp = -f , X G n 

(2.5) V-u = 0, xGfi 

(2.6) u = 0, xedQ. 
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The spatial domain occupied by the fluid is denoted by Q. In the SELM approach [2], the 
forces acting on the elastic structures are accounted for through a force density acting on 
the fluid given by 

(2.7) f(x,f)=AF. 

The F denotes the collection of forces acting on the elastic structures. The A denotes the 
fluid-structure force coupling operator of the SELM approach, see [2]. 

A key property of which we shall make use is the commutation of the Laplacian operator 
A and the divergence operator V-. This can be expressed as 

(2.8) A(V-) = (V-)A. 

By taking the divergence of equation |2.4| and using equation |2.5[ we have 

(2.9) Ap = V ■ f . 
Formally, this has the solution 

(2.10) p=--A-V-f. 
Substituting this for p in equation |2 .4| yields 



(2.11) /iAu = -pf 

(2.12) p = X-(V)A-^(V-). 

The formal solution for u can be expressed as 

(2.13) u = pM"'A-V(-f). 

We have used the property that p is a projection operator and for solutions u of equation |2.4| 
we have pu = u. 

In the SELM approach the dynamics of the elastic structures is given by [2] 

The effective stochastic dynamics given by equation |2.1| for X is obtained by substituting u 
from equation |2 . 1 3| into equation |2 . 14| and using equation |2.7| for f . This derivation uses that 
the dissipative operator is the Laplacian A which is symmetric, we discuss an alternative 
form for the effective stochastic dynamics of elastic structures when the dissipative operator 
is not symmetric in Appendix [A] 

2.1. Discretization of SELM Hydrodynamic Coupling Tensor. To utilize 
the SELM hydrodynamic coupling approach in practice requires numerical methods for the 
approximate computation of the coupling operator i?sELM- One approach is provided by 
discretizing and numerically approximating solutions of the underlying fluid equations. 

To facilitate the development of the stochastic numerical methods, it will be convenient 
to work with numerical discretizations which have a number of special properties. We shall 
find it convenient to work with linear coupling operators which when discretized satisfy the 
following adjoint condition 

(2.15) A = ar^. 

Such a condition is closely related to the requirement that the fluid-structure coupling con- 
serves energy, see [2] [25]. For the numerical approximation of the projection operator p we 
shall find it useful for the approximating discretized operator p to satisfy p — 0^ . For the 
discretized Laplacian we shall require L = . We shall also find it useful for the p to 
commute with the Laplacian pL = Lp. One realization of the SELM approach having these 
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properties is the Stochastic Immersed Boundary Method, see [4] |25] . These conditions are 
imposed in this initial presentation to avoid a number of technical issues and for clarity in 
the exposition. In practice, it is likely many of these conditions can be relaxed. 
The effective equations for the elastic structures can then be expressed as 

(2.16) ^=^SELM(F)+g 

at 

(2.17) HsELM = -VpL-^pK. 

The F denotes now the discretized collection of forces acting on the elastic structures. For 
the discretized model to be consistent with statistical mechanics it is required that 

(2.18) G = (gg^) = 2fcsTi?sELM. 

This covariance structure can be found as a consequence of the principle of detailed balance 
applied to the discretized system, see [^. 

3. Generation of the Stochastic Driving Fields Accounting for Thermal 
Fluctuations. The stochastic driving field g is required to have the covariance structure 

(3.1) G = 2fcflrii'sELM. 

Generating random variates with these statistics is in general made challenging since the 
covariance structure involves correlations which decay like ~ 1/r, for the case of three spatial 
dimensions. Using traditional methods such as Cholesky factorization is expensive. For 
Cholesky factorization the computational cost scales as O(M^), where M is the number of 
rows in G. For systems involving many degrees of freedom, M will be large. 

As an alternative approach, we shall generate g accounting for the thermal fluctuations 
utilizing the underlying discretized fluid equations. The covariance G can be expressed using 
the special form of the discretized hydrodynamic coupling operator Hselm in equation 2.17 
This is given by 

(3.2) G = -KL-^K'^ 

(3.3) K = V2kBTrp. 

A useful consequence of expressing the covariance in this form is that g can be generated 
from a Gaussian random fleld ^ having the covariance structure C = —L~^, in the case L is 
symmetric. We discuss an alternative when L is non-symmetric in Appendix [A] 
When L is symmetric, equation |3.2| provides a method for generating g from 



(3.4) g = Ki. 
The covariance of g is then given by 

(3.5) (gg^) = K{ie)K'^ 

(3.6) (^^^) = C = -L-\ 

This reduces the problem of efficiently generating g to the problem of efficiently generat- 
ing a Gaussian random fleld ^ with the covariance structure G = —L~^. This is challenging 
in general since the covariance structure is the inverse Laplacian, which has row entries which 
decay from the diagonal with a scaling like ~ 1/r, for the case of three spatial dimensions. 
In the notation, for the (i,j)-entry we take r = |i — j|Aa;. For the random field ^, this 
corresponds to long-range spatial correlations which decay like ~ 1/r. 
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3.1. Stochastic Iterative Methods for Generating Correlated Variates. 

To generate the random field ^ with the required long-range spatial correlations, we shall 
exploit the special property of the covariance structure that it is obtained as the inverse of a 
sparse matrix. As in solving linear systems of equations, this property is suggestive that an 
iterative approach may be efficient provided we can generalize such iterative methods to the 
stochastic context. For lattice theories, stochastic iterative methods have been introduced 
which generalize traditional iterative methods such as SOR, Gauss-Siedel, and Jacobi itera- 
tions to generate random fields, see [2 [5^. This was further generalized to obtain stochastic 
multigrid iterative methods in [131 114) . 

For the SELM approach, we shall develop stochastic multigrid iterative methods for 
generating the random variates ^. This will allow for the efficient generation of the stochas- 
tic driving fields g in the effective hydrodynamic equations for elastic structures given in 
equation |2.1| by using equation |3.4| We now discuss this stochastic iterative approach for 
generating random variates in more detail. 

To generate the stochastic driving fields we shall develop a Gibb's sampler for the re- 
quired multi-variate Gaussian distribution. A Gibb's sampler is developed with stochastic 
iterations constructed in a manner which exactly preserves the target distribution as an in- 
variant measure of the iterative process. The efficiency of the Gibb's sampler is determined 
by two important factors. The first is the computational cost required to compute each 
iteration of the sampler. The second is the number of iterations required to obtain a random 
variate having negligible correlations with the previously generated random variates. 

For multi-variate Gaussian distributions, we shall use the following linear stochastic 
iterations 



(3.7) Z"+' =7?,Z''-Hs-^r?". 

The T/'"' are taken to be independent Gaussian random variables with mean zero and co- 
variance J, 

(3.8) (r,") = 

(3.9) {r/'"(r7")^) =5^,„J. 



In the notation, 5m,n denotes the Kronecker (S-function, the superscript T denotes the vector 
transpose, {•) denotes a probability expectation. The stochastic iteration given in equa- 
tion |3.7| can also be expressed in terms of the probability density p"(z) at iteration n. These 
probability densities satisfy 

(3.10) p"+'(z) = y^7r(z,w)p"(w)dw 

(3.11) Tv(z, w) = , exp [(z — Rw — s)^ J^^ (z — _Rw — s)l . 

V27rdetJ L J 

We now discuss how a stochastic iterative method having the form of equation |3 . 7| can 
be used to sample a multi-variate Gaussian with a specified mean fi and covariance C. For 
this purpose, expressions can be derived which relate the iteration matrix R and s to the 
mean /i and covariance C, [131 IT4] . The mean fj, is given by 

(3.12) n^{I-Ry^s. 

This follows from equation |3.7[ by taking the expectation of both sides and taking the limit 
as n — > oo. The covariance of Z"'^^ is given by 

(3.13) C"+' = ((Z"+i - /i)(Z"+' - fif). 
This satisfies the following linear recurrence equation 



(3.14) 



C"+' = RC"R'^ + J. 
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This follows from equation 3.13 by using equation 3.7 to express Z""""^. Taking the limit as 



oo of equation |3. 14| we obtain 

(3.15) J ^ C ~ RCR^. 

This provides a condition on the covariance J for the stochastic driving field which 
is necessary for the iteration based on R to exhibit the target covariance C. A sufficient 
condition for the invariant measure of equation |3.7| to be unique can be obtained by re- 
writing equation |3.15| as 

(3.16) AC = J. 

Retaining the matrix indexing when treating C as a vector, the entries of this operator can 
be expressed as 

(3-1''') -^(»l,»2),(jl .J2) ~ ^il,i2^jl,j2 ~ -^!l.jl-^^2,j2■ 

The invariant measure is then ensured to be unique provided the linear operator A acting 
on C is non-singular (invertible) in equation |3.16[ This provides a condition on R. 

We now discuss for a specified mean /x and covariance C a stochastic iterative scheme 
for sampling the multi-variate Gaussian. For the given covariance structure C, equation |3.15| 
gives the required covariance J for the stochastic driving field. To generate the stochastic 
driving field, we shall find a factor Q so that 

(3.18) J = QQ'^. 
The random variates r] can then be generated using 

(3.19) 'n = Qn 

where n is a vector having components which are each independent Gaussians with mean 
zero and variance one. The computational efficiency of this approach will depend on the 
particular form of Q and whether a sparse factor can be found satisfying equation |3.18| 

To make this procedure concrete, we shall consider a stochastic sampler based on Gauss- 
Siedel iterations. In the deterministic setting, the Gauss-Siedel iteration approximates solu- 
tions of the following linear system 

(3.20) At, = s. 

The A is assumed to be symmetric. In the Gauss-Siedel iteration, the solution is approxi- 
mated by splitting the matrix into a diagonal part D, upper triangular part U , and lower 
triangular part L so that 

(3.21) A = D - L-U. 



For symmetric A we have U = . An iteration of the form of equation 3.7 is then performed 
with 

(3.22) R^{D-Ly^U. 

We now show that a stochastic iterative method can be formulated based on such Gauss- 
Siedel iterations for sampling a Gaussian with a specified mean |x and covariance C = A~^ . 
As in the deterministic case, the computation of Rz can be performed especially efficiently 
provided A = is sparse. However, in the stochastic setting it is also required that Qn be 
computed each iteration, which could be computationally expensive depending on the factor 
Q of equation |3. 18| which is used. 
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We now discuss one such approach for obtaining a factor Q satisfying equation|3.18| For 



Gauss-Siedel iterations, the covariance J of the stochastic driving term in equation 3.7 takes 
the specific form 

(3.23) J = QQ^ ^{D~L~ U)'^ - (D - L)'^U{D - L - U)'^L{D - L)"^. 
An expHcit factor for Q can be found [14] and is given by 

gg^ = iD- L)-^ [(D-L- u){D -L- uy\D - U) 

(3.24) + U{D -L- Uy\D - L - U)] (D - Uy^ 

= {D - Ly'^ D{D - uy^ . 

Since U = and the transpose of the inverse is the inverse of the transpose, we have that 
the following factor Q satisfies equation |3.23| 

(3.25) Q ^ (D - Ly"^ D"^^^ . 

From this explicit form, we see that Q can be computed efficiently provided A is sparse. 
This follows since (D — L)^^ is a lower triangular matrix whose inverse can be found using 
back-substitution. The D is the diagonal matrix whose square root is readily computed. 

In the case that the covariance C has a sparse inverse, we have that A — is sparse. 
When A is sparse with a constant number of non-zero entries per row, Qn can be computed 
with a computational costs of only 0{N), where A*' is the number of rows of A. For C with 
sparse inverse, one stochastic iteration of the sampler has a total computational cost of only 
0{N). 

The second cost associated with the Gibb's sampler is the number of iterations required 
to obtain a new random variate which has a negligible correlation with the previously gener- 
ated random variates. We now compute the autocorrelation of the random variates generated 
by the stochastic iterative method given in equation |3.23[ The correlation between the vari- 
ate generated at iteration and iteration k is given by 

(3.26) C'"'*"' = ((Z° -/x)(Z'' -/x)^). 
This can be shown to satisfy the linear recurrence equation 

(3.27) (^(o,fc) ^^fc^(o,o)^ 



This follows by using equation 3.7 to express in terms of Z" and r]-' [14] , 

This indicates that the number of iterations required for the correlations between random 
variates to become small depends on the initial correlation and the decay rate of the matrix 
R. This provides the following bound in terms of matrix norms 

(3.28) ||C(°''='|| < ||i?||'=||C<'''°'||. 

To make more precise what is meant by the correlation becoming negligible between random 
variates, we shall require the following bound to hold 

(3.29) ||C<''''='|| < 6 

where e is small. The number of iterations k required to obtain random variates which have 
negligible correlation is then bounded below by 

(3.30) fc>log(E/||C(°-«)||)/log(||ii||). 

This analysis shows there is a direct link between the number of iterations required to 
obtain a negligible correlation and the number of iterations required in traditional determin- 
istic iterative methods to obtain a high level of accuracy in approximating the solution of a 
linear system. This suggests that techniques which improve the convergence of traditional 
iterative methods also should improve the sampling efficiency of stochastic iterative methods. 
One widely used approach which dramatically improves the convergence of iterative methods 
is the use of multigrid methods. We explore this approach in detail in the next section. 
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Fig. 3.1. Multigrid operators. In multigrid, three operators are defined for each level of refine- 
ment of the mesh. The first is a smoother operator which is used to iteratively approximate the 
solution of the linear system of equations on the mesh at level I. For transmission of information 
between meshes at different levels of resolution, a restriction operator X^~^ and a prolongation op- 
erator X^_-^ are defined. The restriction operator X^~^ maps data from a more refined mesh at level 
I to data on a less refined mesh at level £ — 1. The prolongation operator maps data from a 

less refined mesh at level £—1 to data on a more refined mesh at level £. These operators are applied 
in combination to obtain iterations in the multigrid method. For more details, see Algorithms 
and\3\ 



3.2. Stochastic Multigrid Methods for Generating Correlated Variates. 

In the deterministic setting, multigrid can be used to greatly accelerate the convergence of 



iterative methods [8] [21]. As can be seen from equation 3.28 reduction in the number of 
iterations required for convergence in the deterministic setting corresponds to a reduction in 
the number of iterations required to obtain random variates which are negligibly correlated. 
We now discuss how the multigrid method can be utilized as a sampler for random variates. 

In multigrid iterations three fundamental operators are utilized, see Figure |3.1[ The 
first is a smoother operator which is used to iteratively approximate the linear system of 
equations on the mesh at level £. For transmission of information between meshes at different 
levels of resolution, a restriction operator and prolongation operator are defined. 

The restriction operator I^"^ maps data from a more refined mesh at level £ to data on a 
less refined mesh at level l—l. The prolongation operator maps data from a less refined 
mesh at level ^ — 1 to data on a more refined mesh at level I. 

When using multigrid for sampling random variates, the primary modification will be 
to use a smoothing operator S which is stochastic. For our present purposes we will take 
the smoothing operation to be the stochastic Gauss-Siedel iteration defined by equation |3.7[ 
equation |3.22| and equation |3.25[ To ensure the target Gaussian with mean /i and covariance 
C is the invariant measure of the stochastic multigrid iterations, we shall also require the 
prolongation operator and restriction operator preserve the variational structure of the linear 
equations ^14j. This variational property is discussed in more detail below. 

For the prolongation operator /|_iin three dimensions we shall use tri-linear interpola- 
tion. For the restriction operator we use the adjoint operator given by 



(3.31) ir = {n-.y 



This ensures that a consistent variational principle is satisfied for the linear systems of 
equations associated with each of the levels of resolution used in the multigrid iterations. 
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The variational principle satisfied at the most refined level of the mesh is given by the solution 
of the linear system of equations being a minimizer of the energy defined by 



(3.32) 



E{v) 



-V Av- 



V b. 



In the multigrid iterations, the smoother at level 
linear system of equations with 



approximates the solution of the 



(3.33) 
(3.34) 
(3.35) 



A'^'^^lfAlf 

b(^) = /rb. 



Provided condition |3.31| h olds, t he solution on the mesh of level £ is the constrained minimizer 



of the energy in equation 



3.32 



over all vectors of the form v 



This follows since 



the energy of equation |3.32 can be expressed as 
(3.36) S(/fvf')) 



A 



)-(-"')' 



We have used equation 



3.31 to obtain the expression for A^^\ 
This variational property ensures that at each level of refinement the stochastic smoother 
samples a Gaussian which is the marginal probability distribution of the target multi-variate 
Gaussian given on the most refined level of the mesh. This follows since for each level of the 
mesh, the smoother samples the probability distribution 



(3.37) 



p(^)(v(^') 



1 



V27rdet4('^) 



exp 



1 



The variational property can be shown to be sufficient to ensure the probability distribution 
of the target multi-variate Gaussian with mean /i and covariance C is the invariant measure 
of the stochastic multigrid iterations [14] . 

As a basis for our stochastic sampler for random variates, we shall use the formulation 
of multigrid referred to as "Fast Adaptive Composite Mesh Multigrid" (abbreviated FAC- 
multigrid) |15l 122) . The detailed steps of the FAC-multigrid are summarized in Algorithms 
[l||3] In the case that j4 is a sparse matrix with only a constant number of non-zero en- 
tries per row the FAC-multigrid iterations can be carried-out with a computational cost of 
0{N \og{N)) operations. In the deterministic setting, it can be shown that FAC-multigrid 
converges with an error with a specified threshold e in 0(1) number of iterations 8, 21; 22j. 
The stochastic sampler based on FAC-multigrid then provides a method to generate random 
variates with negligible correlation with a computational cost of only 0(A'' log(A'')) opera- 
tions. Since FAC-multigrid can be used on domains with boundaries and on spatially adaptive 
meshes, this provides an efficient approach for generating ^ and the stochastic driving term 
g in equation |2.16| 
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Algorithm 1: v ^ FAC-Multigrid(v, b, /xi, /X2)- 
Data: An initial value of v, the right-hand-side b, the number of smoother 

iterations (z/, /ii, /i2). 
Result: An approximate solution v of the linear system = b. 

Procedure: 

1. Compute the residual of the initial value r = — b. 

2. Initialize the initial value for the correction q ^ 0. 

3. Perform a full sweep of the mesh cells on the multilevel mesh 
q Full-Sweep(q, r, v, fii,(J,2)- 

4. Correct the solution v v -|- q. 



Algorithm 2: q'^^ ^ Full-Sweep(q^'^\ r^^), i/, /zi, /Z2). 

Data: An initial value q*^^-* , the right-hand-side r^^' , the number of smoother 

iterations {v, fii, ^12)- 
Result: An approximate solution q^^^ of the linear system A^^W = r^^^. 

Procedure: 

1. If the current mesh is at the level of refinement common to all meshes in the 
hierarchy then perform one V-Cycle on the mesh: 

qW ^ V-Cycle(qW,rW,/ii,/i2). 

2. Otherwise, coarsen the residual to the next mesh level r^^"^' I^~^r^^K 

3. Perform a full sweep of the cells of the multilevel mesh 
q«?-i) ^ Full-Sweep(q(^-i),r(^-i),zy,/ii,/i2). 

4. Apply the correction to the solution on the current level: q'*^' ^ /|_iq^^~^^- 

5. Apply the smoother with the initial value q'^' for v iterations for the linear 
system A^^W = r^^^. 
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Algorithm 3: ^ V-Cyclc(v(^), b^). 

Data: An initial value v^^), the right-hand-side h^^\ the number of smoother 
iterations (/ii,/Z2)- 

Result: An approximate solution v*^^^ of the linear system A^^'^w = h^^\ 
Procedure: 

1. Apply the smoother with the initial value v^^' for ni iterations for the linear 
system ^(''-'w = b'-''-'. 

2. If the current mesh is at the coarsest level of refinement in the hierarchy 
then skip to step 5. 

3. Otherwise, perform a V-Cyclc on the next coarsest mesh in the hierarchy. 
Let b(^-i) ^ (b(^-i) AvW), v('^-i) ^ 0, then compute 

v(«-i) ^ V-Cycle(v(^-i),b(^-i)). 

4. Correct the solution on the current level v^^' <— v'^^ + I^g_i^y^^~^^ ■ 

5. Apply the smoother with the initial value v^^' for ij,2 iterations for the linear 
system A^^W = b*^^). 
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4. Spatially Adaptive Discretizations. For many physical systems the level of 
resolution required to accurately describe the state of the system can vary dramatically over 
the spatial domain. When approximating such systems numerically, a natural approach is to 
use discretizations which have different levels of spatial resolution. Many spatially adaptive 
finite difference discretizations have been developed for this purpose [3 1151 1161 1231 128j . 

A significant challenge arising in the stochastic setting for such discretizations is that the 
discrete Laplacian L is often non-symmetric. The use of symmetry was central in the deriva- 
tion of the effective equations for the elastic structures satisfying the principle of detailed 
balance in Section [2] In this report we present initial results for a widely used discretization 
on adaptive meshes in which the discrete Laplacian L is non-symmetric. 

As an alternative approach to obtaining effective stochastic dynamics for the elastic 
structures, we shall use the stochastic dynamics which would arise for the elastic structures 
from the fluctuations of the underlying discretized fluid equations with discrete Laplacian L 
and with a discrete fluctuation-dissipation principle satisfled. This can be shown formally 
to give the following effective dynamics for the elastic structures (see Appendix [AJ 

(4.1) dlCt ^ HFdt + RdBt 
where 

(4.2) H = -Fpa^i^^i^V^r^ 

(4.3) = -Fp (afj.-^L~^C + CL''^ii~^a^ p'^F^. 

In these expressions the operator L is no longer required to be symmetric. The operator C 
gives the equilibrium covariance structure of the fluctuations of the fluid velocity consistent 
with Gibbs-Boltzmann statistics of the discretized system, see [2] |3] . 

For the equations of hydrodynamics, we shall use a MAC discretization of the Laplacian 



defined at cell centered nodes on a spatially adaptive structured mesh |15II23] . see Figure 4.1 
To obtain such discretizations on multilevel meshes, we express the Laplacian in terms of 
the gradient and divergence operators 

(4.4) A = Laplacian 

(4.5) 2? = V • Divergence 

(4.6) g = V Gradient. 

To approximate the operators, we define for any discretization mesh a partition of the spacial 



domain {fimjm, see Figure 4.1 For a given partition cell Q,^ we allow for numerical values 
to be defined both at the center of the partition cell and at the center of the faces of the 
partition cell. We approximate the Divergence Operator T> at the center of a partition cell 
using 

(4.7) {Dh)^^- ^h^^k - ■arn,k- 

^^■^ fe=i 

The term hm,k denotes the vector value at the center of the fc"' face of the partition cell Q,m. 
The b denotes the composite vector of all such values on the partition. The nm.fe denotes 
the outward normal to the fc*'' face of the partition cell. The term Aa;m is the width of the 
partition cell. The notation {■)m denotes the component corresponding to the value at the 
center of the partition cell with index m. A useful property of this approximation to the 
divergence operator is that its evaluation only requires at the face centers the components 
in the normal direction, see the dot product in equation |4.7| 

We approximate the Gradient Operator Q at the center of the faces of each partition 
cell. Given the different levels of resolution in the mesh, many cases can arise in principle. 
By convention, we restrict our methods to deal with meshes which have the nested property 
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Multilevel Mesh Quad-Tree Data Structure 




Fig. 4.1. Mesh Cells at a Typical Coarse-Refined Interface. 



that neighboring cells differ in resolution by at most one level. This requires only two cases 
be considered at each face of a partition cell. The first is when the neighboring cell is at the 
same level of spatial resolution. This corresponds to Axm = Ax^^, , where £k denotes the 
index of the neighbor in the direction of the fc"* face of the partition cell. The second is when 
the neighboring cells differ by one level of resolution, Axm = lAxi,^, or Asm — \A.xt.f.- 

To approximate the gradient operator on a face shared with a neighbor at the same level 
of resolution, we use 



(4.8) 



(Gc 



= sign(ni, 



In the notation (Om.fc denotes the components corresponding to the vector value at the center 
of the k^^ face of the partition cell with index m. The notation denotes the fe*'' vector 

component. The discrete gradient operator only defines the fc"* vector component at each 



face since this is all that is required by the discrete divergence operator D of equation |4.7[ 

To approximate the gradient operator on faces shared between neighbors differing by 
one level of spatial resolution, we must consider a cluster of partition cells. To simplify the 
discussion, we consider the case where the partition cell with index m has neighbors at the 
fe"* face which are of a more refined level of resolution, Aa;m — 2Aa;£j. . We define the cluster 
to be the collection of partition cells consisting of the partition cell with index m (labeled 
A) and the four neighboring partition c ells in the direction of the outward normal of the 



V'^ face (labeled B,C,D,E), see Figure |4.l[ The components of the gradient operator are 
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approximated by 

(4.9) (Gc)(,^: 

(4.10) (Gc)g=: 

(4.11) (Gc)^^' 

To obtain a discretization of the Laplacian A on meshes with multiple levels of resolution, 
we use the approximation 

(4.12) L = DG. 

The discrete gradient operator G and discrete divergence operator D are defined by equa- 
tions [X7J|4TT] Similar discretizations have been used in |15l 1231 f2H] . 

Using this approach to discretize the Laplacian allows for both Neumann and Dirichlet 
boundary conditions to be imposed readily on rectangular domains. For Neumann conditions 
the domain is discretized so that faces of the partition cells align with the domain boundary. 
To impose the Neumann conditions the values of components of the gradient are specified 
at the center of faces of the partition coinciding with the boundary. For Dirichlet boundary 
conditions the domain is discretized so that the centers of the partition cells align with the 
domain boundary. To impose the Dirichlet boundary conditions the values are specified at 
the center of partition cells coinciding with the boundary. The Laplacian is then computed 
using equation |4.12[ where the range of the gradient and divergence operators are restricted 
to the non-boundary values of the partition cells. 

5. Performance of the Stochastic Samplers in Practice. We now demon- 
strate how the stochastic numerical methods perform in practice for generation of the random 
variates ^. We consider the case of spatially adaptive discretizations, similar results are ex- 
pected for domains with boundaries. The random variates ^ which are sought are required 
to have the covariance structure 

(5.1) C = {^e) = ~2L-'C. 

The L is the linear operator on the spatially adaptive mesh which approximates the Laplacian 
differential operator and is no longer required to be symmetric. It can be shown that for 
the specific spatially adaptive discretization discussed in Section [4] the L~^C is symmetric. 
For the purposes of comparison, the stochastic multigrid approach is compared with the 
stochastic Gauss-Siedel iteration, see Sections |3. 1| and |3.2[ 

As a measure of the efficiency of the generation methods we consider for a given number 
of iterations the strength of the correlation between random variates generated with each of 
the stochastic iterative methods. We consider the correlations of random variates for spatial 
discretizations both in two and three spatial dimensions, see Figure It is found that the 
stochastic multigrid method greatly reduces the correlations in random variates generated by 
the stochastic iterative methods. In fact, using stochastic multigrid appears to only require 
0(1) number of iterations to generate random variates with negligible correlation. The use 
of Gauss-Siedel stochastic iterations along show a high level of correlations even after many 
iterations. This is found both for two and three spatial dimensions. 

These numerical studies show that the stochastic multigrid method provides a very ef- 
ficient approach for generating random variates with long-range correlations on spatially 
adaptive meshes. The presented approach provides a method to generate such random vari- 
ates with a computational complexity of 0{N \og{N)) number of operations. This provides 
a dramatic improvement over Cholesky factorization approaches, which require 0{N^) op- 
erations to generate the factors and 0{N^) to generate each random variate. 



■ sign(n 



(fc) 



I (Cs + Cc) - CA 



■ sign(n 



(fc) I (CB + cc) - CA 



cl'« + (Gc)g'] 
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These initial results indicate that the developed stochastic multigrid methods provide 
a potentially versatile tool for generating the stochastic driving fields in hydrodynamically 
coupled systems. We have demonstrated here an approach by which stochastic multigrid 
methods can be developed which allow for the simulation of fluid-structure interactions 
in the SELM formulation with a computational cost of only 0{N log{N)) operations. The 
developed real-space stochastic numerical methods allow for simulations on spatially adaptive 
meshes and on flow domains having boundaries. The stochastic numerical method for the 
stochastic dynamics of elastic structures with SELM hydrodynamic coupling is summarized 
in Algorithm |4] 
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Fig. 5.2. Comparison of samplers based on Multigrid and Gauss-Siedel. For spatially adaptive 
meshes in two and three spatial dimensions, a comparison is shown when using the multigrid sampler 
and Gauss-Siedel sampler. For the multigrid based samplers nearly independent random variates ^ 
are generated by the stochastic sampler after only a few iterations. For Gauss-Siedel many iterations 
are required to yield nearly independent random variates ^. 
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Algorithm 4: Stochastic Dynamics of Elastic Structures with Hydrodynamic Cou- 
pling 

Data: Parameters for the physical model and numerical methods. 

Result: Stochastic dynamics of the hydrodynamically coupled elastic structures. 

Procedure: 

1. Compute the forces F acting on the elastic structures. 

2. Compute the velocity of the elastic structures by evaluating 

-ffsBLM — ~rpct-^~^pr"^) which uses the discretization of the fluid equations 
i, p and the fluid-structure coupling operator F. 

3. Generate the stochastic driving term g accounting for thermal fluctuations. 
This is done by using the stochastic multigrid method (Algorithm [I]) to 
generate the random field ^ with covariance — L^^C. The stochastic driving 
term is generated by g = Fp^, see equation |3.4| 

4. Update the configuration of the elastic structures using an SDE temporal 
integrator (i.e. Euler-Maruyama Method). 

5. Return to step 1. to compute the next time step for the dynamics of the 
elastic structures. 
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6. Conclusions. This technical report shows a proof-of-concept approach for efR- 

cicntly generating the stochastic driving fields in fluid-structure systems using the SELM 
formalism for hydrodynamic coupling. In future work, a more detailed numerical investi- 
gation will be presented demonstrating the efficiency of the presented methods for spatial 
domains having more complex geometries and for specific discretizations developed for the 
projection operator on spatially adaptive meshes. Many of the basic ideas presented here are 
expected to generalize to be applicable more broadly in the development of stochastic nu- 
merical methods for the efficient generation of stochastic fields with long-range correlations 
often required in the simulation of spatially extended stochastic systems. 
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Appendix A. Derivation of effective stochastic dynamics of elastic struc- 
tures. We now derive formally a set of effective equations for the dynamics of the elastic 
structures. The equations only involve the elastic structure degrees of freedom and eliminate 
the fluid degrees of freedom. For this purpose, we consider increments of the form 

(A.l) AX(f) =X(t + Af)-X(f). 

An increment over the time [t,t -\- At] of the stochastic dynamics for the elastic structures 
X is given by 

/t+At 
Vu{q)dq. 

In the limit in which the fluid relaxes to a quasi-steady-state on the time scale of the elastic 
structure dynamics, it is expected that X will be an Ito or Statonovich stochastic process. 
Formally, the statistics of this process can be determined by the mean and covariance of its 
increments [l2l[24| . 

For this purpose, we denote the mean of an increment of X by 

We denote the covariance of an increment of X by 

(A.4) b = ((AX- AX) (AX- AX)^)/At 

where AX = (AX) . A stochastic process for the effective dynamics of X having increments 
with approximate mean a and covariance b is given by 



(A. 5) dXt = a(t,Xi)dt + fT(t,Xt)dBt. 

where cra^ = b. Here we assume an Ito limiting process, but a Statonovich limiting process 
may also be considered. A more rigorous justification of how to remove the fluid degrees of 
freedom can be obtained by considering a singular perturbation analysis of the Backward 
Kolomogorov Equations, see |20j . 

To determine a and b, we shall derive formally the leading order expressions in terms of 
e and At, (with the assumption that e <^ At <gi 1). To simplify the notation we consider the 
case when t = Q without the loss of generality. Throughout, we shall use the approximation 

/•At rAt 

(A.6) AX= / ru{q)dq^r / u(q)dq. 

Jo Jo 
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The solution to the fluid equations can be expressed as 
(A.7) 

Jo Jo 
For the mean, the leading order term is given by 

rAt _j 

(A.8) a = r / e}'-"''' '^f{q)dq 

Jo 

(A.9) = -Ater£.-^t(0) + o{At). 

This gives 

(A.IO) a(t,X) = -erz:"'f(o). 

We can express the covariance as 

(A.ll) ((AX- AX) (AX- AX)^) =rQ*r^ 

(A. 12) * = / dr ds(l>{r,s) 

Jo Jo 

(A.13) <^(r,s) = {u(s)u^(r)). 

Using equation |A.7| and equation |A.8| the covariance can be expressed as 
(A. 14) v[( = *i + *2 

rAt rAt 



rat rat _j 

(A.15) *i = / dr dae"' ^C(O) 
Jo Jo 



At rAt , „ , . 

e-l£{r-9). 



(A.16) *2 = / ds I dr f [ ^^'''-''>Q{dBgdBt,)Q' e' 
Jo Jo Jo Jo 

TheC(O) — (u(0)u"^(0)) and is independent of At and e. The integral for v^i can be performed 
analytically to obtain 

(A.17) *i = e""/:-^ (^e'^''~''^ - C(0) (e^*^"'^^ - £.-'^ . 

Under the assumption e ^ At this gives a leading order term of order e^. 

To compute "^2 we use the Ito Isometry, which can be expressed formally as (dB^dB^) = 
S{q — w)dqd'w. This yields 



(■At rAt 

(A.18) 4'2 = / ds dr I dw I dq ''''•''-''> QS{q - w)Q' e 



10 JO Jo Jo 

rAt rAt rrAs 

(A.19) = ds dr dw e ''■^''-'^'QQ' e' 

Jo Jo Jo 

(A.20) ^h + h 

(A.21) Ji = ds r dr f dw e'"''^('-'"'OQ^e'"'^'^'"-'") 
Jo Jo Jo 

(A.22) h = ds dr f dw e'^"'^('"-™'gQ V"'^'^'''-'"' . 

Jo Ja Jo 

By changing the order of integration we can evaluate Ji to obtain 

^At ^Af r3 _ _ 

(A.23) /i = / dw ds dr '^^'"'^^ QQ^ 
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(A.24) = Ai + 

(A.25) Ai = ££-1 / dw ds e' ^'"-"'gQ^e' ^ 

(A.26) = eL-^QQ^L-'^ (^e/I"^ (^e-^*'"'^"^ ~ ^) ~ ^*^) ' 

Under the assumption e ^ At we have Ai is of order ^ . 
By changing the order of integration in A\ , we obtain 

(A.27) Ai=e£-Vi 

(A.28) Ji=/ ds / ^(^-"''QQ^e'= ^ 

Next, we use that QQ^ = — e^^£C — C£^e^^, this allows for Ji to be expressed as 
(A.29) Ji = ds dt« A 

(A.30) = Ate - e / dr e'^^Ce"^ 

(A.31) = Ate - ^e£-^ [e^'^^'^Ce^*'"'^^ - c] . 

From equation |A.27[ the leading term when e At is given by 
(A.32) Ai = Ate£-^C(1 + o(l)). 

From equation |A.23| this gives the leading order expression for Ji 
(A.33) h = Ate£"^C(l + o(l)). 

A similar analysis can be carried out to show that 
(A.34) /a = AteC£"^(l + o(l)). 

From equation |A.14| this shows that to leading order 
(A.35) * = Ate(^£"'C + C/:"^) (l + o(l)). 

The covariance b of the increments is then given from equation | A. 11| by 
(A.36) b(X,t) = ear (^r^'C +C£"^) r^. 

When C^^C is symmetric, this simplifies to 

(A.37) b(x,t) = -2ear/:~'cr^. 

This formal analysis suggests using for the effective stochastic dynamics of the elastic 
structures 

(A.38) dXt = HFdt + RdBt 

where 

(A.39) H = -Va^l-^ C'^V'^ 

(A.40) RR^ = (aij.-^C-^C + CZI^^^^'a) F^. 
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In the case that L = pA and C = ksTii^ — y) this gives the same effective stochastic 
dynamics for the elastic structures as equation |2.1| Since pA = Ap and pC = Cp, we can 
express the above operators as 

(A.41) H = -raAt"VA"V^r^ 

(A.42) RR^ = -rp {ayr^^~^C + CA^^^i^'a) p^F^. 

For an approximation of the Laplacian A by the discrete operator L, this provides the 
foUowing formal approximation of the operators appearing in the elastic structure equations 

(A.43) H = -VpoLpT^L'^p^v'^ 

(A.44) RR^ = -rp [onT^L'^C + CL^'^/i'^aj p^F^. 

The utility of this last expression is that it was derived without the assumption that L 
is symmetric, only that L have negative eigenvalues. The formal derivation provides one 
approach for obtaining effective stochastic dynamics for the elastic structures from the un- 
derlying fluid equations when the discretization of the Laplacian is non-symmetric. As we 
discuss in Section [4] this is especially important when using discretizations developed for the 
Laplacian on spatially adaptive meshes which are often non-symmetric. 

We should emphasize the above analysis is quite formal. The derivations presented here 
are meant primarily to motivate the use of the reported stochastic dynamics for the elastic 
structures. These expressions should be treated with some caution and ultimately need to 
be established either through a more rigorous analysis or through numerical validation. 



